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When a classical system is driven through a continuous phase transition, its nonequilibrium 
response is universal and exhibits Kibble-Zurek scaling. We explore this dynamical scaling in the 
novel context of a three-dimensional topological magnet with fractionalized excitations, namely the 
liquid-gas transition of the emergent mobile magnetic nronopoles in dipolar spin ice. Using field¬ 
mixing and finite-size scaling techniques, we place the critical point of the liquid-gas line in the 
three-dimensional Ising universality class. We then demonstrate Kibble-Zurek scaling for sweeps 
of the magnetic field through the critical point. Unusually slow microscopic time scales in spin ice 
offer a unique opportunity to detect this universal nonequilibrium physics in current experimental 
setups. 


I. INTRODUCTION 

The theory of equilibrium phase transitions is one of 
the major achievements of 20th-century physics [1]. Ef¬ 
fective theories based on local order parameters describe 
the universal aspects of such transitions. Our under¬ 
standing of nonequilibrium physics in their vicinity has 
also progressed on many fronts, from scaling properties 
of dynamical correlation functions [1, 2], to Kibble-Zurek 
(KZ) behavior [3-5]. 

Topological phases do not fit into the conventional 
framework as they are not characterized by local order 
parameters [6]. Indeed, one of the principal attractions 
of systems with emergent gauge fields is that standard 
Landau-Ginzburg-Wilson arguments can fail. The phe¬ 
nomenon of deconfined quantum criticality [7] and anal¬ 
ogous phenomena noted earlier in quantum dimer mod¬ 
els [8] are examples of this. On general grounds this sug¬ 
gests that critical points can exhibit new and unexpected 
universality classes, and it is not clear in which settings 
this happens. 

This work combines the study of a critical point in a 
topological system with an investigation of topological 
phases out of equilibrium. This is a young, multi-faceted 
and rapidly developing subject [9]. In particular, the 
interplay of nonequilibrium physics, long-range interac¬ 
tions and classical topological order near a critical point 
is largely unchartered territory, which this work explores. 

We consider three-dimensional dipolar spin ice, 
realized in rare-earth pyrochlore oxides such as 
{Dy,Ho} 2 Ti 207 . At low temperatures, frustration pre¬ 
vents the system from ordering [10] and it enters a 
highly degenerate topological regime described by an 
emergent gauge field with magnetic Coulomb-interacting 
monopoles as fractional excitations [11]. In many— 
primarily thermodynamic—respects, spin ice is well 
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modeled as a magnetic version of an electrolyte [12], 
whereas- particularly out of equilibrium- fundamental 
deviations on account of the Dirac strings connecting the 
monopoles have been observed [13, 14]. 

In an appropriately oriented magnetic field, spin ice 
exhibits a liquid-gas phase diagram in the temperature- 
field plane, where a first-order transition line terminates 
at a critical point [15]. This behavior is extremely un¬ 
usual for a system of localized spins. However, its origin 
is naturally understood in the monopole picture: the field 
acts as a tunable chemical potential for the monopoles, 
which form a type of Coulomb liquid [16], well known to 
exhibit liquid-gas phase diagrams [17]. 

Our study combines the above questions on critical 
and nonequilibrium properties of topological systems by, 
firstly, characterizing this critical point and, secondly, 
studying the universal dynamics in its vicinity. 

We characterize equilibrium properties via field mixing 
analysis [18] and finite-size scaling. This analysis strongly 
suggests that the critical point is in the three-dimensional 
(3D) Ising universality class, as is believed to be the case 
for conventional Coulomb liquids [17, 19]. 

We then analyze local (single spin-flip) dynamics in the 
vicinity of this transition. We study the magnetization 
response on linearly sweeping the system across the criti¬ 
cal point. The response is hysteretic and, for slow enough 
sweeps, universal a la Kibble-Zurek [3, 20, 21]. We estab¬ 
lish the presence of KZ scaling [22-25], confirming Ising 
universality and allowing us to obtain the dynamical scal¬ 
ing exponent. 

Although there is some experimental evidence for KZ 
scaling of defect density [5, 21], a decisive test of the 
scaling of dynamical response functions in this context 
is still lacking. Our simulations suggest that KZ scaling 
in spin ice is realistically accessible in field sweep exper¬ 
iments [26], owing to unusually slow microscopic time 
scales of the large rare-earth moments [27, 28]. In ad¬ 
dition, the (uniform) magnetization turns out to be a 
direct measure of the monopole density near the critical 
point—an informative quantity that is otherwise difficult 
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FIG. 1. (Color online) Liquid-gas phase diagram of spin ice 
in the (T, H) plane. The insets show representative spin con¬ 
figurations of two tetrahedra in each of the phases. The field- 
and temperaturelike directions at the critical point, H and T, 
are shown. Note that the crossover from kagome ice to spin 
ice occurs outside the range of parameters shown. 

to access experimentally. As such, the liquid-gas tran¬ 
sition is a classical instance of destabilizing the emer¬ 
gent magnetic vacuum via the “Schwinger mechanism” 
of monopole-antimonopole pair creation [29]. 

Spin ice may thus be a unique experimental system, 
not only in which to observe this KZ scaling, but also as 
an instance of a 3D topological material out of equilib¬ 
rium. 


II. FIELD MIXING AT THE CRITICAL POINT 


ice, monopole excitations are activated and exponentially 
sparse at low T, whereas in the saturated phase they are 
dense as they correspond to the lowest-energy spin con¬ 
figuration. H acts as a chemical potential and the transi¬ 
tion is characterized by a finite jump in monopole density, 
from a monopole gas to a monopole liquid phase [15, 16]. 
Representative spin configurations of two tetrahedra in 
each of the phases can be seen in the insets of Fig. 1. 

The first-order line terminates at a critical point, typ¬ 
ical of liquid-gas phase diagrams. To determine its loca¬ 
tion, it is convenient to introduce the reduced Hamilto¬ 
nian 

j; = -/?£ - hm, (2) 

where /? and h are dimensionless couplings; m is the 
dimensionless [111] component of the magnetization; 
and e is the dimensionless total spin interaction en¬ 
ergy. The location of the critical point in the (T, 77) 
plane is determined from the crossing of the Binder 
cumulant of the magnetization m for different system 
sizes L (see Appendix B for details). Our best esti¬ 
mate for the critical point thus obtained is (T c , H c ) = 
(0.5875 ± 0.0005 K, 0.86295 ± 0.00005 T). 

In liquid-gas phase diagrams, the fieldlike and temper¬ 
aturelike directions at the critical point typically do not 
coincide with the original parameters in the system (H 
and T here). The field mixing formalism [18, 31], de¬ 
veloped in the context of classical fluids, allows the uni¬ 
versal properties of such critical points to be extracted. 
According to the revised scaling hypothesis, the correct 
scaling operators rh and e are linear superpositions of the 
quantities m and £ that appear in the Hamiltonian, 


We consider dipolar spin ice in a magnetic field of 
strength H along the [111] direction. We use the Hamil¬ 
tonian 


U = -fxH ^2 (H • <Ji + J ^2 

i (ij) 


+ °rLJ2 

j>i 
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where {c^ = ±1} are pseudospins of local easy axis 
e ? ; and magnetic moment [i = 10 J = —1.24 K 
and D = 1.41 K are respectively the nearest-neighbour 
exchange and dipolar coupling constants relevant for 
Dy 2 Ti 2 C >7 [30]; r nn is the nearest-neighbour spacing of 
the pyrochlore lattice; r I? is the separation vector be¬ 
tween sites i and j, and H = FH is the external mag¬ 
netic field which we measure in Tesla (T). We set k-g = 1 
and measure all energies in Kelvin (K). 

The schematic phase diagram of spin ice in the (T, H) 
plane is shown in Fig. 1. At low T, upon increasing the 
field strength H , spin ice crosses over to a partially po¬ 
larized phase known as kagome ice, and then undergoes 
a first-order transition to a saturated phase. In kagome 



where r and s are appropriate mixing parameters [18]. 
The scaling fields conjugate to m and e—H and T 
respectively—define new fieldlike (symmetry-breaking) 
and temperaturelike (non-symmetry-breaking) directions 
in parameter space. The symmetry of the critical point 
and the universal scaling content of the transition is man¬ 
ifest in the mixed operators and fields. 

From the locations of the Binder cumulant minima 
(alternatively, the maxima of the susceptibilites), the 
slope of the first-order line in the vicinity of the criti¬ 
cal point is determined. This slope in the (/3, h) plane, 
(d/3 c /dh)\ h _ h , is directly related to the mixing parame¬ 
ter r = (d/? c /d/i) _1 = 4.283 ± 0.005. It is interesting to 
note that this value is in reasonable agreement with the 
estimate obtained from the Clausius-Clapeyron relation, 
dH c /dT = —AS/AM ~ 4.45 [32], where AM and AS 
are, respectively, the differences in magnetization and en¬ 
tropy between the kagome ice and saturated phases. The 
relation is expected to hold at low T where the transition 
is strongly first order and has indeed already been shown 
to extrapolate up to the critical point in 16. 
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FIG. 2. (Color online) Finite-size scaling collapse of the equi¬ 
librium susceptibility \ i n the vicinity of the critical point, 
using the critical exponents 7 CT and v a as fitting parameters. 
The collapse gives 7 CT = 0.76 ± 0.02 and v a = 0.41 ± 0.01, in 
agreement with the Ising values 7 ^ ~ 0.79 and v a ~ 0.40. 
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FIG. 3. (Color online) Histogram of mixed operator fluctu¬ 
ations Am, Ae, where AX = X — {X ), normalized by their 
respective standard deviations, -J((AX) 2 }. The character¬ 
istic shape is indicative of emergent Z 2 symmetry and Ising 
universality. Inset: Corresponding data for the 3D Ising and 
3D XY models at criticality (adapted from Ref. [33]). 



Obtaining the mixing parameter s is not as straight¬ 
forward [33]. Here we do so by requiring the statistical 
independence of fluctuations in to and e at the critical 
point, and find s = 0.0(3) ± 0.0(7). Other methods give 
consistent results (not shown). 

A substantial admixture of rh in e is expected from the 
monopole picture. Indeed, the spin interaction energy 
in the Hamiltonian encompasses both a Coulombic term 
as well as a chemical potential for the monopoles, and 
they are of comparable strength. While it is possible to 
separate (approximately) the two contributions in e, we 
find that it does not lead to an appreciable improvement 
in the field-mixing analysis. 

Performing a finite-size scaling analysis in the vicinity 
of the critical point we obtain critical exponents that are 
consistent with 3D Ising universality: 7 ^ = 0.76 ± 0.02 
and v a = 0.41 ±0.01, compared with the Ising values 
7 c- ~ 0.79 and v a ~ 0.40 in nonzero field [40]. The finite- 
size scaling collapse for the magnetic susceptibility x is 
illustrated in Fig. 2. In the 3D Ising class, fluctuations in 
to are dominant over those in e and we observe the same 
magnetic susceptibility exponent 7 ^ whether we consider 
fluctuations in to, to, or even e. In order to observe 
the heat capacity exponent substantial statistical 
accuracy would be required to ensure that e does not 
contain any contribution from to. 

The joint probability distribution of the fluctuations 
in to and e at the critical point, the form of which is 
known to be universal [34, 35], is another powerful tool to 
identify the universality class of a system [18, 33, 36, 37]. 
The distribution can be seen in Fig. 3. It displays a 
characteristic shape which is the hallmark of an emergent 
Z 2 symmetry and thus Ising universality. The data can 
be compared with the 3D Ising model (Z 2 symmetry) 
and contrasted with the 3D XY model [[/(l) symmetry] 
in the insets of Fig. 3. Direct comparison between Ising 


and XY criticality in three dimensions is important due 
to the very similar scaling exponents, which are difficult 
to differentiate between within our numerical accuracy. 

Recently, the long-standing question [17] of the uni¬ 
versality class of critical classical Coulomb liquids has 
been resolved in numerical simulations [19] in favour of 
Ising (and not mean-field) behavior. Our results are con¬ 
sistent with this, in accordance with the picture of a 
liquid-gas transition of the emergent monopoles. It is 
important to note that the Ising universality does not 
arise trivially from the Ising nature of the original spins. 
Indeed, their Z 2 symmetry is explicitly broken by the 
applied field. The spins constitute a vacuum for the 
magnetic monopole excitations. It is then this emer¬ 
gent Coulomb liquid, or magnetolyte, which undergoes 
a liquid-gas transition as would a classical electrolyte, 
with associated emergent Coulombic criticality. This is 
highly nontrivial: emergent monopoles are connected by 
a network of Dirac strings. These strings are however 
statistically and energetically immaterial, thus allowing 
the thermodynamic properties of the system to be the 
same as for a gas of real pointlike charges. 


III. KIBBLE-ZUREK SCALING OF 
HYSTERESIS 

The standard KZ choice is to vary the temperaturelike 
parameter identified above (see e.g., the recent study 38). 
However, in the context of a liquid-gas transition, fine 
tuning would be necessary to identify the precise T(T, H) 
trajectory, which is in general rather difficult away from 
the limit \H — H c \, |T — T c | —> 0. As the simplest ex¬ 
perimental prototcol, we therefore propose to vary the 
applied field H(t), which is then a combination of H and 
T, and measure the magnetization m(t). The resulting 
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FIG. 4. (Color online) Hysteresis in the magnetization m dur¬ 
ing a forwards-backwards field cycle. The field H is swept at 
constant temperature from 0.75 T, through the critical point, 
to 0.95 T in time tq (measured in MC steps per spin), return¬ 
ing similarly. The system size is L = 10. Inset: Hysteresis 
loop size Amhyst(t) as a function of time to the critical point 
t. 

behavior is controlled by the most relevant correlation 
length in the vicinity of the critical point (the magnetic 
correlation length in the case of 3D Ising universality, as 
discussed below). 

Let H(t) be varied linearly from Hi to Hf in a time 
tq such that H( 0) = H c . Here, and in the following, we 
use single spin-flip dynamics and measure time in Monte 
Carlo (MC) steps per spin. Let the instantaneous re¬ 
laxation time of the system be £t(i). Initially, £ 4 (i) is 
short and the evolution is adiabatic; the magnetization 
is thus close to its equilibrium value at H(t). However, 
the relaxation time diverges near the critical point, as 
does its rate of change. When its rate of change be¬ 
comes larger than that of the system parameter (of order 
one for linear sweeps), the system can no longer stay in 
equilibrium: (d£t/di)|t=t KZ ~ 1. This identifies the KZ 
time, Ikz- Crudely speaking, for t £ [—t kz^kz]) the sys¬ 
tem is out of equilibrium, the dynamics are slow, and 
the magnetization response lags behind the field. For 
t fi<Z; the evolution is once again adiabatic. Over a 
full (forwards-backwards) sweep cycle, one thereby ob¬ 
tains a hysteresis loop in the magnetization m as a func¬ 
tion of H. Examples of such hysteresis loops for different 
ramp times tq are shown in Fig. 4. The faster the sweep 
(the smaller tq), the larger the hysteresis loop. 

For slow enough ramps (tq >-1), generalized KZ scal¬ 
ing relations predict that the nonequilibrium contribu¬ 
tion to the magnetization assumes a universal form [25]: 

<m(i)> ~ (^) ' tKZ = T ^’ (4) 

where A is the scaling dimension, and v = (3 — A) -1 the 
correlation length exponent, of the most relevant oper¬ 
ator with projection on m close to the critical point; z 



FIG. 5. (Color online) Kibble-Zurek scaling collapse of 
the hysteresis loop size Amh ya t(f), using the critical expo¬ 
nents v and z as fitting parameters. The collapse gives 
v = v a = 0.42 ± 0.01 and z = 1.85 ± 0.05, consistent with the 
Ising values. Both the L — 8 and 10 data are plotted. 

is the dynamical exponent, and Q is a universal scaling 
function. 

The size of the hysteresis loop—namely, the difference 
between forwards and backwards curves, Amh ys t(t)—can 
be calculated as a function of time to the critical point t 
(inset of Fig. 4). Figure 5 shows the KZ scaling collapse 
of Amhyst(t) according to Eq. (4), using the exponents v 
and z as fitting parameters. We find v = = 0.42 ±0.01 

and 2 = 1.85 ± 0.05, again consistent with the values 
for the 3D Ising model: z ~ 2.0 for Metropolis single 
spin-flip dynamics [39]; and (in nonzero field) the most 
relevant operator is the magnetization, with correlation 
length exponent v a ~ 0.40 [40]. Similar results are ob¬ 
tained for linear sweeps in generic directions in the (T, H) 
plane. We also checked robustness to small variations in 
the choice of the location of the critical point. 

In order to observe universal KZ scaling in numeri¬ 
cal simulations of finite-size systems, one has to be in 
the appropriate speed regime [38]. The system can gen¬ 
erally be characterized by three length scales: a lattice 
scale a, the KZ length scale Ikz ~ t K £, and the system 
size L. Universal KZ physics characteristic of the ther¬ 
modynamic limit appears when a <C Ikz L. Too slow 
sweeps lead to Ikz ^ L and adiabatic evolution at finite 
size. For the larger systems considered here (L > 8), 
the universal scaling behavior is clearly visible and spans 
approximately two decades in tq (see Appendix C). 

It was recently demonstrated [14] that spin ice systems, 
as magnetic Coulomb liquids, exhibit a (transient) second 
Wien effect, whereby the monopole density increases in 
response to an applied field. While this is in principle 
relevant to the dynamics discussed above, the effect is 
appreciable only at low monopole densities, which is not 
the case in our work. 

Finally, it is interesting to note that, even during 
our fastest sweeps, the spins located between adjacent 
kagome planes remain essentially fully polarized so that 
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the monopole density p follows extremely closely the 
magnetization m of the system [41] according to the 
formula m(p) oc (5 p+ 10)/3 (in units of pb per spin). 
The magnetization is therefore an excellent proxy for the 
monopole density [42]. 

IV. OUTLOOK AND EXPERIMENTS 
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We have presented a study of the universal equilib¬ 
rium and out-of-equilibrium properties of an emergent 
liquid-gas critical point in a three-dimensional topolog¬ 
ical magnet—spin ice in a [111] field. This holds the 
promise of experimental verification. One ingredient of 
practical importance is an unusually slow microscopic 
time scale due to the large energy barriers of the single¬ 
ion crystal-field environment, combined with comparably 
small transverse fields [43]. Spins appear to flip at a char¬ 
acteristic rate of approximately 1 kHz [27, 28]. This is 
to be contrasted with typical magnetic materials where 
microscopic dynamics occur on time scales of the order 
of nanoseconds or even picoseconds. 

Moreover, single spin-flip Monte Carlo dynamics have 
proved to capture reasonably well real dynamics in spin 
ice materials [28, 44]. Combining these two observa¬ 
tions to translate our results into experimentally rele¬ 
vant terms, we find that KZ scaling in spin ice may be 
accessible by sweeping fields in the range of 0.5-1 T, at 
temperatures of the order of 0.6 K, at rates from around 
0.7 to 0.03 T/s, while measuring magnetization with an 
accuracy of about 1 % or < 0.1 pb per spin. This is emi¬ 
nently accessible and indeed comparable to achievements 
in earlier field sweep measurements [26]. 

Spin ice may thus be the ideal experimental setup to 
observe dynamical Kibble-Zurek scaling in the context of 
a three dimensional topological magnet with a liquid-gas 
critical point driven by emergent fractionalized monopole 
excitations. 
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Appendix A: Technical details of the Monte Carlo 
simulations 

Our Monte Carlo simulations of spin ice use a con¬ 
ventional cubic unit cell containing 16 spins. A system 
consists of L x L x L unit cells, so N = 16L 3 spins in 
total, with periodic boundary conditions. The Hamilto¬ 
nian of the system is given by Eq. (1) in the main text 
and we treat the long-range dipolar interactions using 
Ewald summation [45, 46]. We use the experimentally- 
determined parameters relevant for Dy 2 Ti 207 , namely 
J = —1.24 K, and D = 1.41 K from 30, and we imple¬ 
ment Metropolis single spin-flip updates. 

In our out-of-equilibrium simulations, we assume that 
Monte Carlo steps represent the actual time evolution of 
the system, namely, that the system has a well-defined 
single spin-flip time scale. This has been argued to cor¬ 
rectly capture the dynamics in experiments at the tem¬ 
peratures relevant to this work [28]. 


Appendix B: Location of the critical point 


As discussed in the main text, we analyze the fourth- 
order (Binder) cumulant [47] to determine the location 
of the critical point. This cumulant is defined for any 
operator X as 


Y = ((AX) 4 ) 

4 “ ((AX) 2 ) 2 ’ 


(Bl) 


where AX = X — (X). The Binder cumulant of the 
magnetization B™, at constant temperature T, shows a 
minimum at some value of the field which is the finite- 
size value for the transition field at that temperature and 
system size L. Fig. 6 shows the minimum values of B™ as 
a function of the temperature T for different system sizes 
L up to L = 12 (27 648 spins). The crossing for different 
system sizes is indicative of a continuous transition and 
the temperature at which it occurs gives an estimate for 
the critical temperature T c . We find the crossing point of 
the two largest pairs of system sizes to be the same within 
errors (suggesting that finite-size effects for these system 
sizes are not the dominant source of error), and equal to 
T c = 0.5875 ± 0.0005. The estimate for the critical field 
H c is the transition field at T c as indicated by either the 
susceptibility maximum or Binder cumulant minimum. 
Our estimate for the location of the critical point is thus 
(T C ,H C ) = (0.5875 ± 0.0005 K, 0.86295 ± 0.00005 T). We 
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FIG. 6. (Color online) Crossing of the fourth-order (Binder) 
cumulant of the magnetization, B™(T), for different system 
sizes L. The behavior is indicative of a continuous transition 
at T c = 0.5875 ± 0.0005 IC. 


remark that the value of the Binder cumulant close to 
the critical point, B™ = 1.86 ± 0.02, differs from the 
infinite-volume Ising value exactly at criticality, B™ ~ 
1.60. Deviations in the value of the Binder cumulant 
have been reported in the literature due to details of the 
lattice structure and interactions [48, 49]. 

We note that the magnetization m used to obtain the 
Binder cumulant is not the critical Ising magnetization¬ 
like operator rh due to field mixing (discussed in the main 
text). However, identifying m requires the values of T c 
and H c , which in turn would in principle need the Binder 
cumulant of rh. Our analysis is justified because fluctu¬ 
ations of rh dominate close to the critical point, so m or 
any quantity containing some rh can be used. For con¬ 
sistency, we repeated the Binder cumulant analysis using 
rh after obtaining the mixing parameters and checked 
that both the location of the critical point, and the value 
of the Binder cumulant at criticality, remain unchanged 
within error bars (not shown). 


Appendix C: Further analysis of the ramp speed 
regimes 

As discussed in the main text, in order to observe 
the universal out-of-equilibrium behavior characteristic 
of the thermodynamic limit in our sweeps, the correla¬ 
tion length when the system falls out of equilibrium, Zkz, 
should be much larger than the lattice constant a (so 
that the dynamics are universal) but significantly smaller 
than the system size L (so that the evolution is not adi¬ 


abatic) [25, 50-52], 

To investigate the different regimes systematically, we 
measure the total nonequilibrium contribution to the 
magnetization over a whole forwards-backwards sweep 
cycle, i.e., the integrated hysteresis loop area. Using the 
finite-time scaling relation given by Eq. (4) in the main 
text, the total area of the hysteresis loop scales with ramp 
time tq as 

/ Am h yst ( t ) dt ~ T q 2 ^ ^ 1,2+1 ^. (Cl) 

Fig. 7 shows the total hysteresis loop area as a function 
of ramp time tq, for different system sizes L. Departures 
from power-law scaling can be seen at small and large tq . 
The departure at large tq is because the evolution is adi¬ 
abatic when Zkz is comparable to the system size. Indeed, 
this departure occurs at smaller values of tq for smaller 
L, as predicted by Ikz ~ L. The departure at small tq is 
because Zkz ~ a and the evolution is nonuniversal. Con¬ 
sistently, this departure is independent of system size. 
The nonuniversality in this fast sweep speed regime also 
manifests itself in an appreciable direction dependence 
in the nonequilibrium contribution to the magnetization 
(not shown). The universal scaling window is in the in¬ 
termediate tq regime. This window is appreciable (ap¬ 
proximately two decades) for the largest system sizes we 
can access. The pink dotted line in Fig. 7 corresponds 
to scaling with the critical exponents from Fig. 5 in the 
main text. 



T Q 

FIG. 7. (Color online) Hysteresis loop area as a function of 
ramp time tq. Deviations from power-law scaling are seen at 
large tq ~ L z due to the time-evolution becoming adiabatic, 
and small tq due to lattice-scale effects independent of L. 
For intermediate tq there is a universal scaling regime (pink 
dotted line with exponents from Fig. 5 in the main text). 
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